#!/usr/bin/python
"""

t_gene_coding_structure.py
    t_locus
    t_transcript
    t_gene

"""

import os,sys
# add self to search path for testing
if __name__ == '__main__':
    exe_path = os.path.split(os.path.abspath(sys.argv[0]))[0]
    sys.path.append(os.path.abspath(os.path.join(exe_path,"..", "python_modules")))
    myname = os.path.split(sys.argv[0])[1]
    myname = os.path.splitext(myname)[0];
else:
    myname = __name__

import general_util
from general_util import die_error
from wrap_text import wrap_text
from reverse_complement_dna import reverse_complement
from translate_by_codon import translate
from adjacent_pairs_iterate import adjacent_pairs_iterate
from custom_log import custom_log

import copy


status_log = None




from sequence_locus import *
        
        
        
        
        

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   t_transcript

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
class t_transcript (object):
    #
    #    N.B. terminal stops are removed from cdna / peptide sequence
    #       (but obviously not from exon structure)
    #
    UTR_CDS_WRONG_ORDER = 0
    EXONS_OVERLAPPING   = 1
    NON_CODING          = 2
    STOPS               = 3
    UNRESCUED_STOPS     = 4 # no refseq to rescue
    TOLERATED_STOPS     = 5 # > 90% identity with refseq
    CODING              = 6
    coding_status_str   = ["non-coding", "stops", "unrescued stops", "matching refseq stops", "coding"]

    #_____________________________________________
    # 
    #   __init__
    #_____________________________________________
    def __init__ (self, cdna_refseq_id, gene):
        self.gene           = gene
        self.rna_loci = list()
        self.utr_loci = list()
        self.cds_loci = list()
        self.utr5_loci = list() 
        self.utr3_loci = list() 
        self.cdna_refseq_id  = cdna_refseq_id
        self.prot_refseq_id  = None
        self.prot_gi         = None
        self.prot_refseq     = None
        self.cdna_gi         = None
        self.cdna_refseq     = None
        self.cdna_accession  = ""
        self.evidence        = None
        self.genomic_cds     = None
        self.cds_translation = None
        self.coding_status   = t_transcript.NON_CODING
        self.refseq_genomic_match = False
        
        
    #_____________________________________________
    # 
    #   get_acc_id_and_sequence
    #_____________________________________________
    def get_acc_id_and_sequence (self, gene_id, acc_id_and_sequence, cdna_id_prefix, 
                                    prot_id_prefix, gene_id_suffix, use_gi):
        # 
        # only take cds without stops
        #
        if self.cds_loci:
            if use_gi:
                acc_id_and_sequence.append([gene_id, 
                                            self.cdna_accession, 
                                            cdna_id_prefix + self.cdna_gi + gene_id_suffix,
                                            prot_id_prefix + self.prot_gi + gene_id_suffix, 
                                            self.genomic_cds, 
                                            self.cds_translation])
            else:
                acc_id_and_sequence.append([gene_id, 
                                            self.cdna_accession, 
                                            cdna_id_prefix + self.cdna_refseq_id + gene_id_suffix,
                                            prot_id_prefix + self.prot_refseq_id + gene_id_suffix, 
                                            self.genomic_cds, 
                                            self.cds_translation])
                                        

    #_____________________________________________
    # 
    #   __cmp__
    #_____________________________________________
    def __cmp__ (self, other):
        if other == None:
            return -1
        if other.cdna_refseq_id == None:
            return -1
        if self.cdna_refseq_id == None:
            return 1
        return cmp(self.cdna_refseq_id, other.cdna_refseq_id)
    #_____________________________________________
    # 
    #   __str__
    #_____________________________________________
    indent = "      "
    def __str__ (self):
        """
            output as string (for debug purposes) Each field on a new line
        """        
        results = ["    cdna_refseq_id=>" + str(self.cdna_refseq_id)]
        # fields
        for field in ("prot_refseq_id", "prot_gi", "cdna_gi", "cdna_accession", "evidence"):
            results.append(t_transcript.indent + "%-15s = %s" % (field, self.__dict__[field]))

        # loci data
        for field in ("rna_loci", "utr5_loci", "cds_loci", "utr3_loci"):
            if len(self.__dict__[field]):
                results.append(t_transcript.indent + "%-15s =" % field )
                results.extend(map(str,  self.__dict__[field]))

        # loci data
        for field in ("cdna_refseq", "genomic_cds", "prot_refseq", "cds_translation"):
            if self.__dict__[field] and len(self.__dict__[field]):
                results.append(t_transcript.indent + ">" + field )
                results.append(t_transcript.indent + self.__dict__[field])

        results.append (t_transcript.indent + "coding_status = %s" %
                                (t_transcript.coding_status_str[self.coding_status]))
                
        # genomic/transcript correspondence
        if self.refseq_genomic_match:
            results.append(t_transcript.indent + "refseq_genomic_match=True")
        else:
            results.append(t_transcript.indent + "refseq_genomic_match=False")

        return "\n".join(results)
        
        
    #_____________________________________________
    # 
    #   do_exons_overlap
    #_____________________________________________
    def do_exons_overlap(self, all_exons):
        """
        returns true if exons overlaps
        assumes list sorted
        """
        for cnt, (i, j) in enumerate(adjacent_pairs_iterate(all_exons, not all_exons[0].strand)):
            if j.beg < i.beg + i.len:
                #sys.stderr.write("overlapping index = %d\n" % cnt)
                return True
        return False

    #_____________________________________________
    # 
    #   get_exon_strings_for_warning
    #_____________________________________________
    def get_exon_strings_for_warning (self, all_exons, utr_set):
        """
            returns list of exon coordinates
            helper function for sort_cds_utr errors
        """
        exon_strings = list()
        for e in all_exons:
            exon_string = str(e)
            if e in utr_set:
                exon_string += " utr"
            exon_strings.append(exon_string)
            
        return "\n". join(exon_strings)

        
    #_____________________________________________
    # 
    #   sort_cds_utr
    #_____________________________________________
    def sort_cds_utr(self):
        if not len(self.cds_loci):
            return
        strand =  self.cds_loci[0].strand
        utr_set = set(self.utr_loci)
        all_exons = self.cds_loci[:]
        all_exons.extend(self.utr_loci)
        all_exons.sort(key = lambda x: x.beg, reverse = not strand)
        
        if self.do_exons_overlap(all_exons):
            self.coding_status   = t_transcript.EXONS_OVERLAPPING
            # print out warning
            exon_strings = self.get_exon_strings_for_warning (all_exons, utr_set)
            #   put "CDS" in utr5_loci
            #   put "UTR" in utr3_loci
            self.utr5_loci = self.cds_loci
            self.cds_loci = []
            self.utr3_loci = self.utr_loci
            status_log.warning("exon loci overlap "
                                "\nexons=\n%s\n" % exon_strings + 
                                str(self.gene))
            return
            

        orig_cds_loci = self.cds_loci
        self.cds_loci = []
        if len(self.utr5_loci) or len(self.utr3_loci):
            status_log.die_error(str(self))
        five_prime = True
        for exon in all_exons:
            if exon in utr_set:
                if five_prime:
                    self.utr5_loci.append(exon)
                else:
                    self.utr3_loci.append(exon)

            # coding (no more 5' prime utrs. Must be cds or 3'utr)
            else:
                five_prime = False
                self.cds_loci.append(exon)
                if not len(self.utr3_loci):
                    continue

                # 
                #   Error: 3' utr comes after cds
                #       utr and cds interspersed!
                # 
                #
                exon_strings = self.get_exon_strings_for_warning (all_exons, utr_set)
                                                                                                                                                                               # 
                self.coding_status   = t_transcript.UTR_CDS_WRONG_ORDER
                #   put "CDS" in utr5_loci
                #   put "UTR" in utr3_loci
                self.utr5_loci = self.cds_loci
                self.cds_loci = []
                self.utr3_loci = self.utr_loci
                
                status_log.warning("utr and cds interspersed "
                                "\nexons=\n%s\n" % exon_strings + 
                                str(self.gene))
                return
                
                
        begs = [all_exons[0].beg, all_exons[-1].beg]
        ends = [all_exons[0].beg + all_exons[0].len, all_exons[-1].beg + all_exons[-1].len]
        self.coding_locus = t_locus(all_exons[0].chromosome, strand, 0,0)
        self.coding_locus.beg = min(begs)
        self.coding_locus.len = max(ends) - self.coding_locus.beg
        
        
        
        
        
        
    #_____________________________________________
    # 
    #  finalise_data
    #_____________________________________________
    def finalise_data (self):
        """
            handle missing protein and transcript ids
            sort out cds and utr

        """
        if self.prot_refseq_id == '-':
            self.prot_refseq_id = None
        if self.cdna_refseq_id == '-':
            self.cdna_refseq_id = None
        
        
    #_____________________________________________
    # 
    #  check if genomic cds has stops
    # 
    #       if it matches the refseq protein sequence and
    #       that has stops then we shall tolerate stops as well
    #_____________________________________________
    def check_coding_status (self):
        """
        
            check if non_coding (no cds_loci)
            or coding (no stops)
            or stops
            or tolerated stops (matching stops in refseq)
            
            Note:
                This is orthogonal to whether the gene is a coding gene
                I.e. ncbi might have annotated it as pseudo
        """
            
        #
        #   non_coding 
        #
        if not len(self.cds_loci):
            # coding status can be NON_CODING, UTR_CDS_WRONG_ORDER, EXONS_OVERLAPPING
            # leave on current status
            if (self.coding_status > t_transcript.NON_CODING):
                die_error("coding status has been set from somewhere else")
            return
            
        if self.cdna_refseq and self.cdna_refseq.find(self.genomic_cds) != -1:
            self.refseq_genomic_match = True

        stop_position = self.cds_translation.find('X')
        if (stop_position == -1):
            stop_position = self.cds_translation.find('*')

        # 
        # no stop 
        # 
        #   remember that terminal stops have already been
        #       removed from sequence
        #
        if (stop_position == -1):
            self.coding_status = t_transcript.CODING
            return 

        self.coding_status = t_transcript.UNRESCUED_STOPS

        # no refseq, can't rescue stops
        if self.prot_refseq != None:

            self.coding_status = t_transcript.STOPS
            #
            #  calculate %id with refseq
            # assumes refseq aligns
            #
            identity = [0,0]
            for i, j in zip(self.prot_refseq, self.cds_translation):
                if j == '*':
                    continue
                identity[i == j] += 1
            identity = identity[1]  *1.0 / (identity[0] +identity[1])
    
            if identity >= 0.9:
                self.coding_status = t_transcript.TOLERATED_STOPS
    
        
        
        
        
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   t_gene

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
class t_gene (object):
    indent = "  "
    def __init__ (self, gene_id, group_label, locus, gene_id_suffix = ""):
        self.cnt_cds_transcripts = 0
        self.gene_name           = None
        self.gene_id             = gene_id
        self.gene_id_suffix      = gene_id_suffix
        self.group_label         = group_label
        self.transcripts         = dict()  # of t_transcript
        self.gene_locus          = locus
        self.pseudo_locus        = None
        self.gene_type           = None
        self.missing_genomic_sequence = False

        self.cnt_missing_prot_refseq= 0 
        self.cnt_mismatched_genomic = 0 
        self.cnt_genomic_stops      = 0 
        self.cnt_known_genomic_stops= 0 
        self.transcript_evidence    = dict()
        for evidence_type in ["matching", "mismatched", "stops", "tolerated_stops"]:
            self.transcript_evidence[evidence_type] = set()
        
        # false if transcripts do not have genomic sequence
        #   cannot match != mismatch
        self.all_mismatched_genomic = False
        self.all_genomic_stops      = False

    #_____________________________________________
    # 
    #  get_coding_loci
    #_____________________________________________
    def get_coding_loci_before_finalise_data (self):
        """

        in preparation for add_genomic_coding_sequence
        do not rely on cnt_cds_transcripts which is set in finalise
        """

        coding_loci = list()
        for t in self.transcripts.itervalues(): 
            if not len(t.cds_loci):
                continue
            coding_loci.extend(t.cds_loci)
        return coding_loci



    #_____________________________________________
    # 
    #  add_genomic_coding_sequence
    #_____________________________________________
    def add_genomic_coding_sequence (self, exonic_coding_sequence,
                                    coding_exonic_sequence_errors):

        # 
        #   sort out exons into the right order (5'->cds->3')
        # 
        for t in self.transcripts.itervalues(): 
            t.sort_cds_utr()

        try:
            for t in self.transcripts.itervalues(): 
                if not len(t.cds_loci):
                    continue

                genomic_cds  = list()

                for locus in t.cds_loci:
                    if not exonic_coding_sequence.has_key(locus):
                        if not coding_exonic_sequence_errors.has_key(locus):
                            sys.stderr.write(str(locus) + "\n")
                            print coding_exonic_sequence_errors
                            print exonic_coding_sequence
                        raise RuntimeError, str(coding_exonic_sequence_errors[locus])
                    exon_seq, checked_chrm = exonic_coding_sequence[locus]
                    if not len(exon_seq):
                        status_log.die_error("No sequence for " + str(locus))
                    if checked_chrm:
                        locus = copy.copy(locus)
                        locus.chromosome = checked_chrm
                    if not locus.strand:
                        exon_seq = reverse_complement(exon_seq)
                    genomic_cds.append(exon_seq)
                # 
                #
                # DEBUGG
                #t.genomic_cds = ":".join(genomic_cds)
                #t.cds_translation = translate(t.genomic_cds.replace(":", ""))

                t.genomic_cds = "".join(genomic_cds)
                t.cds_translation = translate(t.genomic_cds)
                if not len(t.cds_translation):
                    self.dump_and_warn("Empty translated genomic sequence")


                if t.cds_translation[-1] == '*':
                    t.cds_translation = t.cds_translation[:-1]
                    penultimate_codon = int((len(t.genomic_cds)-1)/3)*3
                    t.genomic_cds = t.genomic_cds[:penultimate_codon]
                    


        except RuntimeError, e:
            for t in self.transcripts.itervalues(): 
                if not len(t.cds_loci):
                    continue
                t.genomic_cds = ""
                t.cds_translation = ""
            self.missing_genomic_sequence = True
            self.reason_for_missing_genomic_sequence = e.message            




    #_____________________________________________
    # 
    #   dump_and_die
    #_____________________________________________
    def dump_and_die (self, reason):
        status_log.die_error(reason + "\n" + str(self))

    def dump_and_warn (self, reason):
        status_log.warning("\nWarning:\n" + 
                        wrap_text((reason + "\n" + str(self)), "\t", "") + "\n")


    #_____________________________________________
    # 
    #   __cmp__
    #_____________________________________________
    def __cmp__ (self, other):
        if other == None:
            return -1
        return cmp(self.gene_locus,  other.gene_locus)

    #_____________________________________________
    # 
    #   __str__
    #_____________________________________________
    def __str__ (self):
        """
            output as string (for debug purposes) Each field on a new line
        """
        results = ["gene_id=>" + self.gene_id]
        
        #
        #   main fields
        for field in ("group_label", "gene_name", "gene_locus", "pseudo_locus", 
                        "gene_type", "cnt_cds_transcripts"):
            results.append(t_gene.indent + "%-15s = %s" % (field, str(self.__dict__[field])))

        #
        #   whether sequence maps
        for field in (  "cnt_missing_prot_refseq", 
                        "cnt_mismatched_genomic", "cnt_genomic_stops",
                        "cnt_known_genomic_stops",
                        "all_genomic_stops", "all_mismatched_genomic", 
                        "reason_for_missing_genomic_sequence"):
            if self.__dict__.has_key(field) and self.__dict__[field]:
                results.append(t_gene.indent + "%-15s = %s" % (field, str(self.__dict__[field])))
        #
        #   each transcript
        if len(self.transcripts):
            results.append(t_gene.indent + "%-15s =" % "transcripts")
            results.extend(map(str,  sorted(self.transcripts.values())))
        
            
            
        return "\n".join(results)

    #_____________________________________________
    # 
    #  compare_genomic_cds_refseq
    #_____________________________________________
    def compare_genomic_cds_refseq (self):

        """
        helper function for finalise_data
        """

        if not self.cnt_cds_transcripts:
            return
            

        self.all_mismatched_genomic = True
        self.all_genomic_stops      = True
        for t in self.transcripts.itervalues(): 
            # ignore non coding
            if t.coding_status <= t_transcript.NON_CODING:
                continue

            # 
            #   stops
            #
            if t.coding_status != t_transcript.CODING:
                # tolerated stops
                if t.coding_status == t_transcript.TOLERATED_STOPS:
                    self.transcript_evidence["tolerated_stops"].add(t.evidence)
                    self.cnt_known_genomic_stops += 1
                else:
                    self.cnt_genomic_stops +=1
                    self.transcript_evidence["stops"].add(t.evidence)
                    
            else:
                self.all_genomic_stops = False


            # 
            #   protein refseq
            #
            if t.prot_refseq == None:
                self.cnt_missing_prot_refseq += 1
                # cannot match != mismatch
                self.all_mismatched_genomic = False
                
            else:
                #   
                #   Genomic/transcript correspondence if has refseq
                #
                if t.refseq_genomic_match:
                    self.transcript_evidence["matching"].add(t.evidence)
                    self.all_mismatched_genomic = False
                else:
                    self.cnt_mismatched_genomic +=1
                    self.transcript_evidence["mismatched"].add(t.evidence)

    #_____________________________________________
    # 
    #  check_coding_transcript_ids
    #_____________________________________________
    def set_ids_if_missing (self, t, transcript_id):
        """
            helper function for check_coding_transcript_ids
        """
        if t.prot_gi == None:
            t.prot_gi = transcript_id
        if t.cdna_gi == None:
            t.cdna_gi = transcript_id
        if t.prot_refseq_id == None:
            t.prot_refseq_id = transcript_id
        if t.cdna_refseq_id == None:
            t.cdna_refseq_id = transcript_id
        
    def check_coding_transcript_ids (self):
        """
            helper function for finalise_data
            
            relies on cnt_cds_transcripts
            
            all coding transcripts should have an id 
            if it doesn't have one, we will give it one!
        """
        
        #
        #   assign identifiers to coding transcripts first
        #
        transcript_index = 0
        for t in self.transcripts.itervalues():
            if not len(t.cds_loci):
                continue
            transcript_index += 1
            transcript_id = "gene_%s_transcript_%d" % (self.gene_id, transcript_index)
            self.set_ids_if_missing (t, transcript_id)


        #
        #   then to non-coding
        #
        for t in self.transcripts.itervalues():
            if len(t.cds_loci):
                continue
            transcript_index += 1
            transcript_id = "gene_%s_transcript_%d" % (self.gene_id, transcript_index)
            self.set_ids_if_missing (t, transcript_id)
            
        
        
    #_____________________________________________
    # 
    #  get_pseudos_with_cds
    #_____________________________________________
    def get_pseudos_with_cds (self, pseudo_types):
        """
        if pseudos, add gene_id and gene type
        """
        if not self.cnt_cds_transcripts:
            return
        
        if self.gene_type == "coding":
            return
        
        pseudo_types.append([self.gene_id, self.gene_type])
            
        
        
        

    #_____________________________________________
    # 
    #  set_gene_type
    #_____________________________________________
    def set_gene_type (self):
        """
        Helper function for finalize data
        """

        #   gene_type = "rna", "pseudo sans cds", "pseudo sans ids", "pseudo with cds", 
        #                       "coding with stops", "coding"
        #
        #   default to RNA type
        #

        #
        #   pseudo
        #
        if self.pseudo_locus != None:
            self.gene_type = "pseudo (ncbi identified sans cds)"
            if not self.cnt_cds_transcripts:
                return

            if self.missing_genomic_sequence:
                self.gene_type = "pseudo (from ncbi, cds, no genomic seq)"
                return

            self.gene_type = "pseudo (from ncbi, cds with stops)"
            for t in self.transcripts.itervalues(): 
                if t.coding_status == t_transcript.CODING:
                    self.gene_type = "pseudo (from ncbi, cds sans stops)"
                    return
                if t.coding_status == t_transcript.TOLERATED_STOPS:
                    self.gene_type = "pseudo (from ncbi, cds stops match refseq)"
            return

        #
        #   sans cds = rna
        #
        elif not self.cnt_cds_transcripts:
            self.gene_type = "rna"
            return
            
            
        if self.missing_genomic_sequence:
            self.gene_type = "pseudo (coding without genomic sequence)"

        else:
            self.gene_type = "pseudo (coding with stops)"
            for t in self.transcripts.itervalues(): 
                if t.coding_status == t_transcript.CODING:
                    self.gene_type = "coding"
                    return
                if t.coding_status == t_transcript.TOLERATED_STOPS:
                    self.gene_type = "pseudo (stops match refseq stops)"




    #_____________________________________________
    # 
    #  finalise_data
    #_____________________________________________
    def finalise_data (self):
        
        #
        #   calculate coding status per transcript 
        #
        for t in self.transcripts.itervalues(): 
            t.finalise_data ()

        #
        #   calculate coding status per transcript 
        #
        for t in self.transcripts.itervalues(): 
            t.check_coding_status ()

        #
        #   compare genomic with refseq
        #   check all coding transcripts have protein and cdna gis
        # 
        #     call after cnt_cds_transcripts 
        #
        for t in self.transcripts.itervalues(): 
            if len(t.cds_loci):
                self.cnt_cds_transcripts +=1

        self.compare_genomic_cds_refseq ()
        self.check_coding_transcript_ids ()

        # 
        #   set gene type
        # 
        self.set_gene_type()
                



    #_____________________________________________
    # 
    #   get_coding_exon_list
    #_____________________________________________
    def get_per_gene_exon_indices (self):
        """
        calculate exon indexes: unique in order for entire gene
        
        exons are sorted in order (reversed for -ve strand) and assign ordinal

        For coding genes, only take transcripts which are coding (no stops)

        """
        #
        #   
        #
        gene_exon_loci = set()
        for t in self.transcripts.itervalues():
            if self.gene_type == "coding" and t.coding_status != t_transcript.CODING:
                continue

            exon_loci = list()
            exon_loci.extend(t.utr5_loci)
            exon_loci.extend(t.cds_loci)
            exon_loci.extend(t.utr3_loci)
            exon_loci = merge_overlapping_loci(exon_loci)
            
            gene_exon_loci.update(exon_loci)

        gene_exon_loci = sorted(list(gene_exon_loci))

        if not len(gene_exon_loci):
            self.dump_and_warn ("No exon loci!")

        if not gene_exon_loci[0].strand:
            gene_exon_loci.reverse()
        exon_loci_to_index = dict()

        for i, exon_locus in enumerate(gene_exon_loci):
            exon_loci_to_index[exon_locus] = i
        return exon_loci_to_index
        
    
    def get_coding_exon_list (self, all_exon_loci, gene_id_prefix, cdna_id_prefix, prot_id_prefix, use_gi):
        """
        return list of loci and exon names shared between exons of this gene
            
        For coding genes, only take transcripts which are coding (no stops)
        
        """

        # 
        #    do nothing if either no cds or missing genomic sequence
        #
        if self.missing_genomic_sequence or not self.cnt_cds_transcripts:
            return
            
        exon_loci_to_index = self.get_per_gene_exon_indices()

        gene_id = gene_id_prefix + self.gene_id

        
        # return array of cdna_id, prot_id, exon_id, chromosome, beg, end, strand
        for t in self.transcripts.itervalues(): 
            
            # For coding genes, only take transcripts which are coding (no stops)
            if self.gene_type == "coding" and t.coding_status != t_transcript.CODING:
                continue
                
            if use_gi:
                seq_cdna_id = cdna_id_prefix + t.cdna_gi + self.gene_id_suffix
                seq_prot_id = prot_id_prefix + t.prot_gi + self.gene_id_suffix
            else:
                seq_cdna_id = cdna_id_prefix + t.cdna_refseq_id + self.gene_id_suffix
                seq_prot_id = prot_id_prefix + t.prot_refseq_id + self.gene_id_suffix

            exon_loci = list()
            exon_loci.extend(t.utr5_loci)
            exon_loci.extend(t.cds_loci)
            exon_loci.extend(t.utr3_loci)
            exon_loci = merge_overlapping_loci(exon_loci)
            
            for locus in exon_loci:
                exon_index = str(exon_loci_to_index[locus])
                all_exon_loci.append([gene_id,
                                      seq_cdna_id,
                                      seq_prot_id,
                                      gene_id + "exon" + exon_index,
                                      locus.tab_delimited ()])



    

    #_____________________________________________
    # 
    #   get_acc_id_and_sequence
    #_____________________________________________
    def get_acc_id_and_sequence (self, acc_id_and_sequence, gene_prefix, prot_id_prefix, cdna_id_prefix, use_gi):
        """
        get accession and sequence

            For coding genes, only take transcripts which are coding (no stops)
        """
        
        # 
        #    do nothing if either no cds or missing genomic sequence
        #
        if self.missing_genomic_sequence or not self.cnt_cds_transcripts:
            return
            
        for t in self.transcripts.itervalues(): 
            # 
            #    For coding genes, only take transcripts which are coding (no stops)
            #
            if self.gene_type == "coding" and t.coding_status != t_transcript.CODING:
                continue
            t.get_acc_id_and_sequence (gene_prefix + self.gene_id, 
                                        acc_id_and_sequence, 
                                        cdna_id_prefix, 
                                        prot_id_prefix, 
                                        self.gene_id_suffix,
                                        use_gi)



